!Calcul de la distribution des vitesses
!Le "fond" pour les temps longs est dû au nombre trop petit de particules. Augmenter N.
!Idem "a2ap-2b.f90" mais avec subruns.

module fonctions
implicit none
integer, parameter :: d=2 !dimension
integer, parameter :: N = 20000000 !nombre initial de particules
double precision, parameter :: p = 0.5d0 !probabilité d'annihilation. Attention: il faut p>0 (prévu: N->0). a_2(p=0)=0
double precision, parameter :: intervalle = 0.01d0 !intervalle log10(N/Nnow) de calcul des a_2
integer, parameter :: Nmin = 2 !nombre de particules minimales d'un run
integer, parameter :: Navg = 40000 !nombre de réalisations où toutes les particules sont annihilées pour moyenne
integer, parameter :: Navgsr = 200 !nombre de SUBRUNS où toutes les particules sont annihilées pour moyenne
double precision, parameter :: alpha = 1.0d0 !coefficient de restitution. Laisser \alpha=1.0
integer, parameter, dimension(1:10) :: valeurs = (/ 1,2,5,10,50,100,1000,10000,100000,1000000 /) !N/Nnow tq histo. calculé
integer, parameter :: histodim = 5000 !nb tot pts: histogramme*2. Ajuster => moyenne mobile
double precision, parameter :: a = -200 !Borne inférieure histogramme
double precision, parameter :: b = +200 !Borne supérieure histogramme
double precision, parameter :: pi = 3.1415926535897932384626433832795028841971693993751058209749445923078164
contains

FUNCTION ran2(idum) !Initialize idum with negative integer. 1.5x +plus lent que ran1. ran2 si +de 100M nb aléatoires
INTEGER idum,IM1,IM2,IMM1,IA1,IA2,IQ1,IQ2,IR1,IR2,NTAB,NDIV
REAL :: ran2,AM,EPS,RNMX
PARAMETER (IM1=2147483563,IM2=2147483399,AM=1./IM1,IMM1=IM1-1,IA1=40014,IA2=40692,IQ1=53668,IQ2=52774,IR1=12211,&
          &IR2=3791,NTAB=32,NDIV=1+IMM1/NTAB,EPS=1.2e-7,RNMX=1.-EPS)
INTEGER :: idum2,j,k,iv(NTAB),iy
SAVE iv,iy,idum2
DATA idum2/123456789/, iv/NTAB*0/, iy/0/
if(idum.le.0)then
	idum=max(-idum,1)
	idum2=idum
	do j=NTAB+8,1,-1
		k=idum/IQ1
		idum=IA1*(idum-k*IQ1)-k*IR1
		if(idum.lt.0)idum=idum+IM1
		if(j.le.NTAB)iv(j)=idum
	enddo
	iy=iv(1)
endif
k=idum/IQ1
idum=IA1*(idum-k*IQ1)-k*IR1
if(idum.lt.0)idum=idum+IM1
k=idum2/IQ2
idum2=IA2*(idum2-k*IQ2)-k*IR2
if(idum2.lt.0)idum2=idum2+IM2
j=1+iy/NDIV
iy=iv(j)-idum2
iv(j)=idum
if(iy.lt.1)iy=iy+IMM1
ran2=min(AM*iy,RNMX)
return
end function ran2

!------------------------------------------------------------------------

subroutine rescale(vitesses,Nnow) !enlève la vitesse du centre de masse
	integer, intent(in) :: Nnow
	double precision, dimension(1:Nnow,1:d), intent(inout) :: vitesses
	double precision, dimension(1:d) :: vcm !vitesse du centre de masse
	integer :: l
	vcm = 0.d0
	do l=1,Nnow
		vcm = vcm + vitesses(l,:)
	end do
	vcm = vcm/dble(Nnow)
	do l=1,Nnow
		vitesses(l,:) = vitesses(l,:) - vcm
	end do
return
end subroutine rescale

!------------------------------------------------------------------------
subroutine histonew(taille,data,histoinout)
	integer, intent(in) :: taille !taille=nb.particules
	double precision, dimension(1:taille,1:d), intent(in) :: data !vitesses
	double precision, dimension(-histodim:histodim), intent(inout) :: histoinout !histogramme
	double precision :: step
	integer :: i,j,numero
	step = (b-a)/(2.d0*dble(histodim))
	do i=1,taille
		do j=1,d
			if (data(i,j) > 0) then !j*h=x => j = int(x/h), et on choisit int plutôt que nint à cause de (i+0.5)*step
				numero=int(data(i,j)/step)
			else
				numero=int(data(i,j)/step) - 1
			end if
			if (abs(numero) <= histodim) histoinout(numero) = histoinout(numero) + 1.d0/(dble(taille)*dble(d)*step)
		end do
	end do	
return
end subroutine histonew

!------------------------------------------------------------------------
subroutine ecriture(data,nom)
	double precision, dimension(-histodim:histodim), intent(in) :: data !histogramme
	character*(*), intent(in) :: nom
	double precision :: step
	integer :: i
	step = (b-a)/(2.d0*dble(histodim))
	open(unit=1,file=nom)
	do i=-histodim,histodim
!		if (data(i) > 0) write(1,*) (dble(i)+0.5d0)*step,data(i) !enregistre la distribution f(c).
		if (data(i) > 0) write(1,*) (dble(i)+0.5d0)*step,data(i)/&
		&((1.d0/dsqrt(pi))*dexp(-((dble(i)+0.5d0)*step)**2)) !enregistre f/M.
	end do
	close(1)
return
end subroutine ecriture

end module fonctions


!----------------------------------------------------------------------------------
!----------------------------------------------------------------------------------


program a2ap
use fonctions
implicit none
integer :: seed,i,j,k,l,p1,p2,Nnow
integer :: Nnowsr
double precision :: beta,norme,prodscal,omega,omegamax,xold,dN,dNnow,v2,mu
double precision :: dNnowsr,omegamaxsr,xoldsr
double precision, dimension(1:N,1:d) :: v
double precision, dimension(:,:), allocatable :: vsr
double precision, dimension(-histodim:histodim,1:size(valeurs)) :: histogramme
double precision, dimension(1:d) :: sigma,vchange
logical :: collision
seed = -1731

dN = dble(N)
histogramme = 0.d0
beta = (1.d0 + alpha)/2.d0

!Navg moyennes
!==========================
avgmain: do k=1,Navg
	collision = .true.
	omegamax = 1.d0 !pas besoin de trouver un meilleur omegamax car les premières iterations l'améliorent
	Nnow = N
	dNnow = dble(Nnow)
	xold = 0.d0
! ***************************************************************	
! *** Initialisation des vitesses entre 1 et -1 de façon isotrope
!	do i=1,N
!		distribisotrope: do
!			do j=1,d
!				v(i,j) = 2.d0*dble(ran2(seed))-1.d0
!			end do
!			if (sum(v(i,:)**2) < 1.d0) exit distribisotrope
!		end do distribisotrope
!	end do
! ***********************************************************
! *** Initialisation des vitesses en loi de Puissance en v^mu
! mu = -1.5d0
!	do i=1,N
!		distribisotrope2: do
!			do j=1,d
!				v(i,j) = 2.d0*dble(ran2(seed))-1.d0
!			end do
!			if (sum(v(i,:)**2) < 1.d0) exit distribisotrope2
!		end do distribisotrope2
!	  do j=1,d
!	    v(i,j) = v(i,j)/(sum(v(i,:)**2)**(mu/(2.d0*(mu+dble(d)))))
!	  end do
!	end do
! ************************************************************************
! *** Initialisation des vitesses selon distrib. gaussienne ("Box Muller")
	do i=1,N
		distribgauss: do j=1,d
			v(i,j) = dsqrt(-2.d0*dlog(dble(ran2(seed))))*dcos(2.d0*pi*dble(ran2(seed)))
		end do distribgauss
	end do
!*************************************************************************
	call rescale(v,N)	
	v2 = sum(v**2)/dN
	call histonew(N,v/dsqrt(v2),histogramme(:,1))
	call ecriture(histogramme(:,1)/dble(k),"a2ap-2-1h.txt")

	!1 run
	!==========================
	boucle_principale: do
		!*** Créer au hasard un vecteur unitaire $\sigma$ de distribution isotrope
		normeiso: do !norme du vecteur doit être inférieure à 1 (rejette les vecteurs en dehors de l'hyper-cercle)
			do i=1,d
				sigma(i) = dble(ran2(seed))-0.5d0
			end do
			norme = sum(sigma**2)
			if (norme < 1.d0) exit normeiso
		end do normeiso
		sigma = sigma/dsqrt(norme)

		!collisions
		!==========================
		differentparticles: do !choisit deux particules différentes au hasard
			p1 = int(ran2(seed)*dNnow + 1.)
			p2 = int(ran2(seed)*dNnow + 1.)
			if (p1 /= p2) exit differentparticles
		end do differentparticles
		prodscal = dot_product(v(p1,:)-v(p2,:),sigma)
		if (prodscal > 0.d0) then !collision seulement si les particules sont dans le même demi-espace
			omega = abs(prodscal) !\omega = ||v_{12} \cdot \sigma ||
			if (omega/omegamax > ran2(seed)) then !collision
				if (p > ran2(seed)) then !Annihilation: refait un tableau de vitesses v "sans trous" mais de dimension Nnow-2
					collision =.true.
					Nnow = Nnow - 2					
					dNnow = dble(Nnow)
					recopie: do !remet aux Nnow premières places les vitesses des particules: enlève les trous
						if ((p1 <= Nnow) .and. (p2 <= Nnow)) then
							v(p1,:) = v(Nnow+1,:)
							v(p2,:) = v(Nnow+2,:)
							exit recopie
						else
							if ((p1 > Nnow) .and. (p2 > Nnow)) then
								exit recopie
							else
								if ((p1 <= Nnow) .and. (p2 > Nnow)) then
									if (p2 == Nnow+1) then
										v(p1,:) = v(Nnow+2,:)
										exit recopie
									else
										v(p1,:) = v(Nnow+1,:)
										exit recopie
									end if
								else
									if (p1 == Nnow+1) then
										v(p2,:) = v(Nnow+2,:)
										exit recopie
									else
										v(p2,:) = v(Nnow+1,:)
										exit recopie
									end if
								end if
							end if
						end if
					end do recopie
					v(Nnow+1,:) = 0.d0
					v(Nnow+2,:) = 0.d0
				else !Collision (inélastique)
					vchange = beta*omega*sigma
					v(p1,:) = v(p1,:) - vchange
					v(p2,:) = v(p2,:) + vchange
				end if
			end if
			if (omega >= omegamax) omegamax = omega*1.1d0 !rescale le omegamax
		end if

		!enregistre les distributions
		!==========================
		if (collision .eqv. .true.) then
			collision = .false.
			do i=2,size(valeurs)
				if ( (nint(dN/dNnow) >= valeurs(i)) .and. (nint(10**xold) < valeurs(i)) ) then
					call rescale(v(1:Nnow,:),Nnow) !important d'enlever VCM avant calcul
					v2 = sum(v(1:Nnow,:)**2)/dNnow
					call histonew(Nnow,v(1:Nnow,:)/dsqrt(v2),histogramme(:,i))
					if (i==2) call ecriture(histogramme(:,i)/dble(k),"a2ap-2-2h.txt")
					if (i==3) call ecriture(histogramme(:,i)/dble(k),"a2ap-2-3h.txt")
					if (i==4) call ecriture(histogramme(:,i)/dble(k),"a2ap-2-4h.txt")
					if (i==4) then !préparation des boucles de subrun
						allocate(vsr(1:Nnow,1:d))
						do l=1,Navgsr
							vsr(1:Nnow,:) = v(1:Nnow,:)
							Nnowsr = Nnow
							dNnowsr = dNnow
							omegamaxsr = omegamax
							xoldsr = xold
										
!==========================================================================================
! DEBUT SUBRUNS: boucle recopiée
!==========================================================================================

							!1 run
							!==========================
							boucle_principalesr: do
								!*** Créer au hasard un vecteur unitaire $\sigma$ de distribution isotrope
								normeisosr: do !norme du vecteur doit être inférieure à 1 (rejette les vecteurs en dehors de l'hyper-cercle)
									do j=1,d
										sigma(j) = dble(ran2(seed))-0.5d0
									end do
									norme = sum(sigma**2)
									if (norme < 1.d0) exit normeisosr
								end do normeisosr
								sigma = sigma/dsqrt(norme)

								!collisions
								!==========================
								differentparticlessr: do !choisit deux particules différentes au hasard
									p1 = int(ran2(seed)*dNnowsr + 1.)
									p2 = int(ran2(seed)*dNnowsr + 1.)
									if (p1 /= p2) exit differentparticlessr
								end do differentparticlessr
								prodscal = dot_product(vsr(p1,:)-vsr(p2,:),sigma)
								if (prodscal > 0.d0) then !collision seulement si les particules sont dans le même demi-espace
									omega = abs(prodscal) !\omega = ||v_{12} \cdot \sigma ||
									if (omega/omegamaxsr > ran2(seed)) then !collision
										if (p > ran2(seed)) then !Annihilation: refait un tableau de vitesses v "sans trous" mais de dimension Nnow-2
											collision =.true.
											Nnowsr = Nnowsr - 2					
											dNnowsr = dble(Nnowsr)
											recopiesr: do !remet aux Nnow premières places les vitesses des particules: enlève les trous
												if ((p1 <= Nnowsr) .and. (p2 <= Nnowsr)) then
													vsr(p1,:) = vsr(Nnowsr+1,:)
													vsr(p2,:) = vsr(Nnowsr+2,:)
													exit recopiesr
												else
													if ((p1 > Nnowsr) .and. (p2 > Nnowsr)) then
														exit recopiesr
													else
														if ((p1 <= Nnowsr) .and. (p2 > Nnowsr)) then
															if (p2 == Nnowsr+1) then
																vsr(p1,:) = vsr(Nnowsr+2,:)
																exit recopiesr
															else
																vsr(p1,:) = vsr(Nnowsr+1,:)
																exit recopiesr
															end if
														else
															if (p1 == Nnowsr+1) then
																vsr(p2,:) = vsr(Nnowsr+2,:)
																exit recopiesr
															else
																vsr(p2,:) = vsr(Nnowsr+1,:)
																exit recopiesr
															end if
														end if
													end if
												end if
											end do recopiesr
											vsr(Nnowsr+1,:) = 0.d0
											vsr(Nnowsr+2,:) = 0.d0
										else !Collision (inélastique)
											vchange = beta*omega*sigma
											vsr(p1,:) = vsr(p1,:) - vchange
											vsr(p2,:) = vsr(p2,:) + vchange
										end if
									end if
									if (omega >= omegamaxsr) omegamaxsr = omega*1.1d0 !rescale le omegamax
								end if

								!enregistre les distributions
								!==========================
								if (collision .eqv. .true.) then
									collision = .false.
									do j=5,size(valeurs)
										if ( (nint(dN/dNnowsr) >= valeurs(j)) .and. (nint(10**xoldsr) < valeurs(j)) ) then
											call rescale(vsr(1:Nnowsr,:),Nnowsr) !important d'enlever VCM avant calcul
											v2 = sum(vsr(1:Nnowsr,:)**2)/dNnowsr
											call histonew(Nnowsr,vsr(1:Nnowsr,:)/dsqrt(v2),histogramme(:,j))
											if (j==5) call ecriture(histogramme(:,j)/dble(Navgsr*(k-1)+l),"a2ap-2-5h.txt")
											if (j==6) call ecriture(histogramme(:,j)/dble(Navgsr*(k-1)+l),"a2ap-2-6h.txt")
											if (j==7) call ecriture(histogramme(:,j)/dble(Navgsr*(k-1)+l),"a2ap-2-7h.txt")
											if (j==8) call ecriture(histogramme(:,j)/dble(Navgsr*(k-1)+l),"a2ap-2-8h.txt")
											if (j==9) call ecriture(histogramme(:,j)/dble(Navgsr*(k-1)+l),"a2ap-2-9h.txt")
											if (j==10) call ecriture(histogramme(:,j)/dble(Navgsr*(k-1)+l),"a2ap-2-10h.txt")
											if (j==11) call ecriture(histogramme(:,j)/dble(Navgsr*(k-1)+l),"a2ap-2-11h.txt")
											if (j==12) call ecriture(histogramme(:,j)/dble(Navgsr*(k-1)+l),"a2ap-2-12h.txt")
											if (Nnowsr <= Nmin) exit boucle_principalesr
										end if
									end do
									xoldsr = dlog10(dN/dNnowsr)
								end if

								if (Nnowsr <= Nmin) exit boucle_principalesr
							end do boucle_principalesr
							open(unit=1,file="state")
								write(1,'(" Run:",i5,"/",i5,", SubRun:",i5,"/",i5,", #TotRuns:",i7,", #p:",i4)') k,Navg,l,Navgsr,Navgsr*(k-1)+l,Nnowsr
							close(1)								

!==========================================================================================
! FIN SUBRUNS
!==========================================================================================

						end do !de la boucle sur l des subrubs
						Nnow = Nnowsr
						dNnow = dNnowsr
						deallocate(vsr)
						xold = 9.d0 !pour que le if de enregistre les valeurs ne soit plus vérifié dans ce Run
					end if !du if i==6, i.e. du programme de subrun. Plus rien sur le subrun après ça.
				end if !du if de enregistre les valeurs
			end do !du i=1,size(valeurs)
			xold = dlog10(dN/dNnow)
		end if

		if (Nnow <= Nmin) exit boucle_principale
	end do boucle_principale
	
	open(unit=1,file="state")
		write(1,'(" Run:",i5,"/",i5,", #p:",i4)') k+1,Navg,Nnow
	close(1)
end do avgmain
write(*,*)
write(*,*) "Programme termine."
read(*,*)
end program a2ap